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Abstract 

We generalize a simple Monte Carlo (MC) model for dilute gases to consider the transport 
behavior of positrons and electrons in Percus-Yevick model liquids under highly non-equilibrium 
conditions, accounting rigorously for coherent scattering processes. The procedure extends an 
existing technique [Wojcik and Tachiya, Chem. Phys. Lett. 363, 3-4 (1992)], using the static 
structure factor to account for the altered anisotropy of coherent scattering in structured material. 
We identify the effects of the approximation used in the original method, and develop a modified 
method that does not require that approximation. We also present an enhanced MC technique 
that has been designed to improve the accuracy and flexibility of simulations in spatially-varying 
electric fields. All of the results are found to be in excellent agreement with an independent multi- 
term Boltzmann equation solution, providing benchmarks for future transport models in liquids 
and structured systems. 
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I. INTRODUCTION 


The precise behavior of electrons and positrons traveling throngh matter is of vital im¬ 
portance in many new and established technologies. Applications snch as solar cells j2], 
radiation dosimetry [3|, material pore-size classification and positron emission tomogra¬ 
phy m all reqnire an nnderstanding of the fnndamental physical processes involved, inclnding 
accnrate knowledge of energy deposition, macroscopic behaviors, and loss rates. 

Althongh the behavior of high-energy particles can be simnlated qnite accnrately with 
condensed history techniqnes [6], at low energies it is important to model the individnal 
interactions of particles colliding with the backgronnd material, and thereby monitor discrete 
energy losses and processes that change the nnmber of particles in the system. In the systems 
that we are investigating, the nnmber density of the charged particles is low enongh that 
the Debye wavelength greatly exceeds the dimensions of the the system, which is known as 
the “swarm” limit of an ionized gas j7]. 

One snccessfnl approach to modeling snch systems is by solving the Boltzmann eqnation 
[8], which is an eqnation of continnity in phase space. Often this approach is limited to ide¬ 
alized geometries, dne to complexities in the nnmerical solntion and application of bonndary 
conditions. However, many real-world systems are too complex for snch an approach to be 
effective, and in any case, alternative methods shonld ideally be nsed for verification. 

The pre-eminent alternative is to nse Monte Carlo simnlations, which have been widely 
employed for snch pnrposes IMS] ever since compnters have been powerfnl enongh to im¬ 
plement them [13]. Monte Carlo simnlations are very flexible, and can easily inclnde featnres 
from systems that are qnite difficnlt to model in any other manner, snch as interfacial effects, 
secondary particles, and inhomogenons media. 

Accnrate simnlations of condensed systems mnst inclnde the effects of coherent scatter¬ 
ing, where the incoming electrons and positrons interact with many particles of the system 
at once. This can occnr when the de Broglie wavelength of the low energy electrons and 
positrons is longer than the mean distance between molecnles of the condensed matter na. 
It is common to ignore these effects in Monte Carlo simnlations of liqnids |16j . becanse they 
are nsnally insignificant for electrons and positrons with energies of greater than ~ 10 — 20 
eV. However, to accnrately treat particle transport at low energies, we mnst inclnde these 
collective effects, nsnally by way of the medinm’s dynamic strnctnre factor S (Ak, Au) |T5] . 
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which contains information about the medium’s characteristic allowed transfers of momen¬ 
tum hAk and energy hjjj. The dynamic structure is exactly what is measured by coherent 
neutron scattering experiments such as d- 


The present study will describe a new Monte Carlo implementation that models struc¬ 
tured matter using a static structure factor, S (Ak), which is an integrated form of the 
dynamic structure factor. We first use a technique described by Wojcik and Tachiya |T] 
to incorporate the static structure factor into our Monte Carlo model, but we assert that 
this method is only accurate for a certain subset of cases. We subsequently extend this 
technique to overcome its limitations. The Percus-Yevick static structure factor is a simple 
analytic static structure factor [18] that can be used as a benchmark to verify the accuracy 
of our simulation. We have performed simulations of a number of Percus-Yevick systems at 
a range of reduced electric field strengths, and we compare our results with those obtained 
by solving the Boltzmann equation detailed in [8]. 

We begin this study with a discussion of the Boltzmann equation approach to coherent 
scattering. We follow this with a brief description of typical Monte Carlo collision mechanics 
for elastic processes in section IIA We then use the Boltzmann equation coherent scattering 


rates to derive a set of modified cross sections in section IIB and identify the approximation 


made in Wojcik and Tachiya’s original method [Tj. In section IIC, we describe a new method 
that we have developed for performing simulations in spatially-varying electric fields. The 


model system that we are studying is extensively described in section IIP and finally, we 
present our results in section |IV[ including comparisons with the results from both our 
implementation of Wojcik and Tachiya’s method, as well as an independent Boltzmann 
equation solution. 


II. THEORY 

A. Coherent scattering 

Designing simulations of swarm transport in liquids and dense gases presents additional 
challenges compared to the ideal gas case. Because the inter-particle spacing of the neutral 
particles is often less than the de Broglie wavelength of the swarm particles, the swarm par¬ 
ticles must often interact with several neutral particles at the same time, which means that 
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any spatial or temporal correlations between said particles will have an effect on scattering 
events. The Cohen-Lekner theory of electron transport na describes these effects in terms 
of two rates of transfer - momentum and energy - that occur independently. 

Cohen and Lekner express the electron distribution function in a basis of spherical har¬ 
monics. They then modify the standard Boltzmann collision integral to include the dy¬ 
namic structure factor S(Ak,a;), as motivated by van Hove’s definition of the ensemble 
cross section [20], and then show that when the necessary integrals have been performed, 
the dependence is only upon the static structure factor. 

Upon solving the equations for the time evolution of the distribution function, they 
ascribe a physical meaning to two of the mean free path lengths that appear in the collision 
integral expansion. The first fully determines the rate of energy transferred from the swarm 
particles. It is independent of the structure of the medium, and is given by the mean free 
path corresponding to single-particle elastic scattering: 


Ao = (noCTm) = no 27 r / dx sin y (1 - cos x) (^sp (e, x) 


( 1 ) 


where no is the number density of the neutral molecules, (Jsp (e, x) is the angle-differential 
elastic cross section for scattering with a single particle (also known as the binary or gas- 
phase cross section), and is the usual definition of the momentum transfer cross section 
in the absence of coherent effects. Throughout the present work, e refers to the relative 
energy in the centre-of-mass frame during a collision, and y represents the angle through 
which the relative velocity is changed. 

The second mean free path partly includes the effect of the medium and contains all 
information about the rate at which momentum is transferred; 


Ai = (nodm) ^ = ^no27r J dx sin y (1 - cos y) asp (e, y) S (Ak)^ 


( 2 ) 


where S (Ak) is the static structure factor as a function of the momentum transferred and 
dm represents a structure modification of the momentum transfer cross section. 

In a recent paper ESI, the explicit rates of energy and momentum transfer were calculated 
with the inclusion of structure in the Boltzmann equation. The components of this transfer 
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due to the collision term, in the case of zero temperature, are: 
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where {) represents averaging over velocity space and n is the number density of the charged 
particles. 

Note that these representative mean free paths should be considered independently, and 
should be only thought of as an average rate of transfer of the relevant quantity, rather than 
as a prescription for separate collision events. We dehne the ratio r(e) = Aq/Ai = dm/cxm- 
In the dilute gas case, r(e) = 1, because the static structure factor of a dilute gas is unity 
for all momentum transfers. However, in a structured medium such as a dense gas or a 
liquid, the ratio can deviate markedly from unity. If r(e) < 1, there is noticeably less 
momentum transfer than in the single-particle scattering case, which can be interpreted as 
a preference towards forward scattering events. In the opposite case of r(e) > 1, more 
momentum transfer occurs, which causes the particle to change direction without losing as 
much energy as it would in the single-particle scattering. 

In the case of an isotropic single-particle elastic cross section, agp (e,x) = (e); as in 


the model described in section III, and the ratio T fe) reduces to: 


T{e) = — = -J dxsmx(l - cosy) 5 ' 



( 5 ) 


where we have assumed the static structure factor depends only on the magnitude of Ak, 
and |hAk| = hAk ~ 2\/2m€sm | in the limit of a small mass ratio m/mo. Throughout the 
present work, m refers to the mass of each charged particle, and mo the mass of each neutral 
molecule. This form of T (e) is sometimes called the angle-integrated static structure factor, 
S (e), and it is this form of the structure factor that is used in several previous works mg. 

In the case of dilute gases, where S (Ak) = 1, the energy and momentum transfer rates 
converge, yielding the single-scattering model in which every energy transfer is accompanied 
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by a momentum transfer. When the transfer rates differ, however, this theory is not directly 
applicable to Monte Carlo modeling because it does not give a microscopic description of 
how much energy and momentum is transferred in each collision between swarm particles 
and neutral particles. 


B. Sampling coherent scattering in Monte Carlo simulations 

Our Monte Carlo simulations are built around sampling sets of scattering cross sections, a, 
that dehne the probabilities of all interactions between the charged particles and the medium. 
Each cross section represents a single type of scattering process, for example elastic, direct 
ionization, or a particular electronic excitation of the neutral. They usually depend upon 
the relative speed of the charged particle during a collision, and on the scattering angle y. 
In the case of a cold background medium, the collision frequency is simply given by: 


V = n„v<TM (v) , 


( 6 ) 


where v is the speed of the charged particle and atot {v) is the integrated sum of all differential 
cross sections at that speed | 22 ]- This quantity is used to stochastically sample the time 


between each collision (further described in section IIC). When a collision is simulated, 
a specihc cross section is randomly selected according to the relative probabilities of the 
available cross sections ESI. In the case of single-scattering collisions with independent gas 
molecules, the amount of energy and momentum transferred is fully determined by the initial 
energy and the scattering angles. 

For structured materials, an approximate theory has been developed by Wojcik and 
Tachiya jTj, who propose a mechanistic model of electron transport in rare gas liquids. In 
what follows, we have extended this model to be more generally applicable to other systems, 
highlighting the approximations and associated errors of Wojcik and Tachiya’s model. 

The presence of structure requires the introduction of additional microscopic processes 
that, at a macroscopic level, produce the same rate of energy and of momentum transfer 
as in the Boltzmann equation formalism detailed in section IIA We choose to do this 
by separating the original, single-particle elastic cross section into three different processes 
depending on the ratio r(e), as illustrated in Fig. These processes have cross sections, 
labeled by which quantities are affected in the collision; Uboth, '^momentum and cXenergy The 


6 




start 



Figure 1. Flowchart detailing how electric fields and coherent scattering are implemented in the 
SSMC code. 


result of a collision from process Uboth is identical to that of a regular single-particle scattering 
collision. For cXenergy; we start with a regular single-particle scattering collision, but set the 
post-collision direction of motion for the particle to be unchanged. This has the effect of 
transferring a minimal amount of momentum whilst maintaining the same energy transfer 
as in (Jboth- For cTmomentum, we perform a regular single-particle scattering collision, but scale 
the post-collision particle speed to be equal to that before the collision. This results in 
exactly zero transfer of energy, but some change of vector momentum. 

The path lengths Aq and Ai in section [II A| correspond to transfer rates of Um = vnoam = 
u/Aq and Um = vnoam = v/Ki for energy and momentum respectively, where v is the speed 
of the charged particle. To achieve these rates, we combine the cross sections in various 
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Figure 2. Schematic diagram of the various elastic cross-sections used in simulating a Percus-Yevick 
liquid {(j) = 0.4). All quantities are given relative to the elastic cross-section for a single particle. 
Note that the crtot > fsp- 

ratios depending on the value of r(e) = Aq/Ai. If r(e) < 1 we wish to decrease the rate 
of momentum transfer, while maintaining energy transfer, and so choose = r(e)(Tsp, 

'^energy = ~ r(e))(Tsp and = 0. Ill the opposite case, r(e) > 1, we achieve an 

increased rate of momentum transfer, by setting = cTsp, cr^omentum = (^(e) — l)crsp and 
'^energy = 0- This gives a total elastic cross section of cxtot = niax(l, r(e))crsp. The complete 
Monte Carlo procedure, which we refer to as the “Static Structure Monte Carlo” (SSMC) 
method, is shown as a flowchart in Fig. 

The procedure outlined above is designed to reproduce the rates of energy and momentum 
transfer in equations Q and (|^. However, it is not obvious that our construction of the 
microscopic processes achieves this goal. In Appendix we show that such a sampling 
process involving these cross sections does indeed satisfy these requirements, to within the 
order of the mass ratio m/mo as mentioned earlier. These differences are small enough 
that they are unlikely to effect electron-atom simulations, though they may be signihcant in 
systems where ions serve as the charged particles. 
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Wojcik and Tachiya [T] studied liquid argon according to the method described above, but 
their mechanistic model, hereafter referred to as the WT method, effectively capped the value 
of r (e) such that it never exceeded unity. This meant that their total collision frequency 
was unaltered from the single-particle scattering case, and simulation of the particles in the 
energy regions where T (e) exceeded 1 could only be considered approximately accurate. The 
difference between the SSMC and WT methods jTj is shown in Fig. where the regions 
labeled (^momentum are absent in their model, and the total cross section modified accordingly, 
so that it is simply agp- For the aforementioned study of liquid argon, such modifications 
were only required in a small energy range for the structure factor that they employed. One 
of the purposes of the present study is to determine how this approximation affects the 
results for a benchmark Percus-Yevick model, where the approximation is more significant. 

C. Precise treatment of electric fields in Monte Carlo simulations 

Electric fields present a particular challenge for this style of Monte Carlo simulation. 
As the collision frequency z/ of a given charged particle is dependent on its energy e (see 
equation]^, the time between collisions r can be altered by the change of energy of the 
particle due to the electric field, even as it is undergoing the transport between collisions. 
Mathematically, the probability of a time between collisions greater than r can be expressed 
as [HI 

P (t) = exp u {e (t)) d?j , (7) 

where the charged particle’s energy is time dependent due to the particle’s passage through 
an electric field. Explicitly performing this integral for every collision would be very com¬ 
putationally expensive, given that the changes in e will depend on the velocity at which the 
particle is traveling and, for non-uniform electric fields, the position of the charged particle 
as well. 

One popular approach uses the method of “null collisions” ini. where the collision fre¬ 
quency is calculated based on the maximum collision frequency z/q that the particle is likely 
to be able to reach during its transport. Using this constant collision frequency, equation ([^ 
can be solved by equating P (r) with a uniformly distributed random number R, in which 
case 

T = —UQ^lnR. ( 8 ) 
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When the collision occurs, a second random number is generated which is used to account 
for this overestiniation by allowing the charged particle to undergo “null” collisions, where 
no exchange of energy or momentum occurs. This procedure suffers from a requirement 
to “backtrack” if the assumed collision frequency is too low, where it must then make a 
second assumption with a higher collision frequency. It is therefore important to minimize 
the number of null collisions and backtracks to optimize the simulation speed, and more 
modern simulations j23] have been designed with this in mind. The null collision method 
effectively amounts to a form of rejection method for sampling from P{t), which means 
that potentially many random numbers are generated for each valid collision. 

The simulation presented here uses an alternative approach. The cross-sections are spec¬ 
ified as a function of energy, but are assumed to be constant within energy bins of width Se. 
These energy bins can be made arbitrarily small, so there is no loss of accuracy provided that 
we are careful to test that the results are independent of the bin width. However, this means 
that it is sufficient to recalculate r only when the energy of the particle changes from one 
bin to another, so until this occurs, the collision frequency in equation (|^ remains constant. 
We therefore design the simulation so that particles can undergo two types of interactions. 
Collisions with neutral particles are described above, and are governed by i, the energy bin 
of the particle. The second type of interaction occurs when the energy bin is judged to have 
changed due to the effect of the electric field. In such an interaction, the only parameter 
that changes is i, which either increases or decreases by one, triggering a recalculation of 
the time until the next neutral particle interaction. The time until the change, t, is deter¬ 
mined by analysis of the particle’s current velocity and the (constant) acceleration that it is 
experiencing due to the electric field. This is given by the smallest real, positive solution to 
the following equation for the kinetic energy: 

(vo -f {at)f = ei± ^6e. (9) 

Recalculating the time until collision r does not require the use of another random number. 
Recalling equation ([^, we now have additional terms for each change in energy bin: 

T = to + ti + ...+ tn 

[In {R) - Uoto - - . . . - Vn-l tn-l] , (10) 

where each is determined by the time required for the particle to change energy from 
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one bin to the next, and each z/j represents the corresponding collision frequency for each 
energy bin at that time. This equation reduces to equation (|^ in the case of constant or 
zero ti- In practice, the simulation maintains a running measure of the remaining “collision 
probability” for each particle. This is the dimensionless quantity in square brackets in 


equation 10 that is divided by the current u to calculate the time until next collision. 


III. BENCHMARK SYSTEM FOR MODELING COHERENT SCATTERING 


A. Percus-Yevick Hard Sphere Model 


To demonstrate and benchmark the new simulation procedure and code, we apply it to 
a simple model system that requires a correct treatment of structured media. One such 
model, frequently used in the literature, is that of a structure for hard-sphere potentials 
obtained by applying the Percus-Yevick approximation as a closure to the Ornstein-Zernike 
equation, which yields a pair-correlation function |25l 126) . which in turn can be transformed 
into a static structure factor via a Fourier transform and directly used in our simulation. 
We use the Verlet and Weiss jT8] structure factor, which includes some corrections to better 
emulate the structure of a real liquid; 


S{Ak) = 1 + 


24:7] 
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( 11 ) 


where t] = cj) — a = f3 = ^ C = ^- This includes a packing density 

parameter, 0, which specihes how closely the hard spheres are packed. It can be written 
in terms of the hard sphere radius r and the neutral number density Uq as 0 = Ivrr^no. 
This structure factor depends only on the magnitude of the momentum exchange during a 
collision. 

We have modeled systems with a range of densities, from </) ^ 0, which approximates a 
dilute gas, to 0 = 0.4, which states that 40% of the volume is excluded by the hard-sphere 
potentials of the neutral molecules. The angle-integrated forms of each of these structure 
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Figure 3. Angle integrated Percus-Yevick structure factors, from equations 0 and ([^, as a 
function of particle energy and volume fraction. 

factors, as described in equation (|^, are shown in Fig. 

B. System parameters 

Our Monte Carlo codes use an event-by-event model, where every collision is considered 
independently. This allows for a great deal of flexibility in the specihcation of scattering 
mechanics, without any of the approximations used by “condensed history” simulations (see 
e.g. |27]). In addition, the swarm approximation — that all transport particles are inde¬ 
pendent — allows the simulation to be run in parallel. This makes it ideal for scheduled 
multi-processor batch jobs, where execution is not necessarily simultaneous or even sequen¬ 
tial. 

We have calculated a number of transport coefficients for comparison with other models. 
The meaning and derivation of all of these coefficients are described in [8] and |2Hj, but a 
short summary is given here. All coefficients are measured as a discrete function of time- 
step ti- During simulation, if a particle’s history crosses a time-step, its properties (eg: 
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position, velocity, position squared) at that time are sampled. We choose the z-axis to be 
aligned with the electric field, since it is the only element of the system that has a preferred 
direction. Each property is added to a separate running total for each time step, and after 
the simulation is complete, the totals are divided by the (in general) time-dependent total 
of the number of particles. This results in the average of a property over all particles, as 
a function of time. While the transport coefficients are in general time-varying, for the 
systems considered in this study all transport coefficients eventually reach a hydrodynamic 
equilibrium after sufficient time has passed. Afterwards, they merely fluctuate in a small 
statistical range about the reported equilibrium value. 

Depending on the required transport coefficients, different properties of the charged par¬ 
ticles must be recorded. In the following table, each definition |28] is of the named property 
at one point in time: 


Transport coefficient 

Definition 

Mean energy 

e = (6) 

Bulk drift velocity 

tr = 1 (r.) 

Bulk longitudinal diffusion 

= i ((r^) - (r^f) 

Bulk transverse diffusion 

Dt = i E...,. i m) - inf ) 


Angle brackets denote an average over all particles, r* represents the appropriate Cartesian 
coordinate of the position of the particle and e represents the energy of the particle. 

In principle, there are two types of transport coefficients, known as “flux” and “bulk”, 
which approximately correspond to per-particle averages and system averages respectively 
|28| . However, because our model contains no non-conservative collision processes, the “flux” 
and “bulk” quantities should be identical, provided enough time samples are taken. We have 
chosen to measure “bulk” quantities where possible, because this implictly averages over 
changes in velocity between time steps, whereas measuring the flux drift velocity and diffu¬ 
sion requires sampling instantaneous values at discrete time values. For the same number 
of time-steps, without any explicit time-averaging of velocity, the bulk quantities have far 
less statistical error. 

Two ranges of reduced field strengths were used. An approximately logarithmic spacing 
of electric field strengths from 0.001 to 100 Td provides a broad picture of the resulting 
behaviors, while a linear spacing of field strengths from 2 to 12 Td provides detail in the 
range in which we expect the two Monte Carlo methods to disagree most strongly. 
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We employ the cold-gas limit, in which the neutrals are considered to be at rest. The 
present method does not accurately support non-zero neutral temperatures, as the static 
structure factor does not contain information about the temperature of the system. We 
are presently formulating a rigorous treatment of non-zero neutral temperatures using the 
dynamic structure factor |29j . 

All of our simulations were performed at different neutral densities as prescribed by the 
volume fraction 0 and the hard-sphere radius of the single particle cross section agp- The 
transport properties presented are independent of the neutral number density, as it scales 
inversely with the electric field strength. In all cases, the hard-sphere cross section for single¬ 
particle scattering was set as asp (e) = 6 A^, while the charged particles were assigned a mass 
equal to that of an electron, m = rUe and the mass of the neutrals was set to mo = 4n. 


IV. RESULTS AND DISCUSSION 

A. Transport coefficients calculated with the new Monte Carlo method. 

The Percus-Yevick hard-sphere system has been previously studied in |2], and a compar¬ 
ison with those results provides a test of our simulation |30]. Figure shows the various 
transport coefficients simulated by the SSMC simulations, using the approach which over¬ 
comes Wojcik and Tachiya’s approximation. It is instructive to compare this figure with 
Fig. to see the field strengths that are affected most strongly by the features of the struc¬ 
ture factors employed. A table of some of our results is in Appendix where we have 
compared them with Boltzmann equation results. The full dataset is available as supple¬ 
mentary material. We have simulated enough particles that in general the Monte Carlo 
statistical error ED is not visible at these scales, being less than 1% in all cases. Agreement 
with the Boltzmann equation results is to within 1% in all cases, so the datasets would not 
be seen as distinct if shown in the above figures, but a detailed comparison is presented in 
section IIV Bl 

The features of the results are discussed in detail in [8j, and we will not repeat that 
discussion in depth here. One key feature is the presence of structure-induced negative- 
differential conductivity (defined in [32]) apparent in the drift velocity: at moderate field 
strengths of 1-10 Td, the drift velocity is inversely proportional to the field strength. This 
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Figure 4. Mean energy e, drift velocity W, and diffusion coefficients and Dt for Percus-Yevick 
model simulations, as a function of reduced electric field E/uq and Percus-Yevick packing ratio </>. 
Error bars are not visible at this scale. 
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Figure 5. Ratio of diffusion coefficients as a function of reduced electric field strength EjuQ and 
packing ratio (j). Note that the ratio for cj) ^ 0 is 0.5 as required. 


is because at low field strengths, the presence of coherent scattering causes an anisotropy in 
particle scattering which allows the particles to be affected more consistently by the held, 
raising their velocity in comparison to the structure-free case. At higher held strengths, the 
mean particle energy is higher, leading to a reduced de Broglie wavelength, which means 
that the charged particles interact with fewer neutral molecules, so the coherent ehects are 
reduced. This results in a net reduction of forward motion despite a higher average energy. 

We would also like to highlight the variation in the anisotropic dihusion as a function of 
0. In the case of a hard-sphere gas with no structure, we expect the ratio Dl/Dt = 0.5 
[7], which, as shown in Fig. is demonstrated by our simulations. When structure is 
introduced, this ratio changes signihcantly. This effect has been previously explored in jSj 
and through the extended Generalized Einstein Relation. We note that a multiterm 
Boltzmann equation solution is required to achieve the accuracy of the SSMC technique. 
The SSMC simulations have no difficulties in accurately representing this anisotropy in the 
velocity distribution function. 
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B. Comparisons with Boltzmann equation and Wojcik and Tachiya’s method. 


In Fig. 1^ we present the percentage difference between our SSMC and WT results 
and the Boltzmann equation solution. Our implementation of the WT method shows good 
agreement over most regions of field strength, however for some larger field strengths, errors 
of up to 5% in the mean energy and up to 35% in the diffusion coefficients become apparent. 
These differences occur when the energies of the electrons are within the regions that are 
truncated by the WT method. We can identify two competing factors that have an effect 
when the particle’s energy is in these regions, causing differences between the SSMC and 
WT methods. Firstly, the collision frequency is enhanced, and since every collision has a 
chance of both removing some energy from the particle and changing the direction of the 
particle away from the direction of the electric field, this means that particles will tend to lose 
energy at a greater rate. However, this is balanced by the presence of the new momentum- 
only collision, which can occur up to 28% of the time in these regions. The presence of such 
collisions will tend to decrease the energy transfer rate. Nevertheless, Fig. [^clearly shows 
that the combination of the effects is observable for the Percus-Yevick 0 = 0.4 case, with a 
peak difference at approximately 8 Td. This corresponds to a mean energy of about 5 eV 
(see Fig. |^, which is at the peak of the F (e) function. That peak is where we would expect 
the WT approximation to be least suitable. 

For the present model, the disagreements with the Boltzmann equation results are less 
than 1% over all field strengths considered. Such differences are of the order of the numerical 
schemes used in the SSMC and Boltzmann equation methods. We suspect the remaining 
differences are a result of energy meshes used in the Monte Carlo codes or Boltzmann 
equation numerical solutions. 


V. CONCLUDING REMARKS 

We have presented a new Monte Carlo simulation code which accurately accounts for 
the effects of structure in non-gaseous systems by employing a modified mechanistic per- 
collision interpretation of the Cohen and Lekner method for solving the Boltzmann equation. 
The SSMC results accurately replicate those calculated via a multi-term solution of the 
Boltzmann equation to within 1%, with agreement significantly improved over those results 
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Figure 6. Mean energy e and longitudinal diffusion Dl percentage difference for each Monte Carlo 
model versus the Boltzmann equation (BE) model, for the Percus-Yevick structure factor at (/> = 0.4. 
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obtained using the WT method, where errors of up to 35% were observed in some transport 
coefficients. Future work will be focused on utilizing a dynamic structure factor to account 
for other structural and collective effects, including non-zero background temperatures. 
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Appendix A: Comparison between Monte-Carlo and Boltzmann equation averages 


As discussed in section IIB the three different processes in the SSMC have been chosen 
in order to reproduce the mean rate of transfer of energy and momentum that was obtained 
in the derivation of the Boltzmann equation. These different processes, corresponding to the 
“cross sections” defined in the main text, are: Uboth, a normal collision that exchanges both 
energy and momentum; cXenergy, a collision without direction change that aims to allow only 
energy exchange; and ^momentum; a collision with only momentum exchange. In this appendix 
we show that, for isotropic scattering from stationary neutrals (i.e. we consider only T = 0 
as in accordance with the regime of validity of the Monte-Carlo simulations), these processes 
give rise to the same rates of energy and momentum transfer as the Boltzmann equation, 
after neglecting terms dependent on the mass-ratio m/mo. Note that in this appendix, 
dashes refer to post-collision quantities. 

In a collision, it is appropriate to consider the particle in the center of mass frame. In this 
frame, isotropic scattering implies that the particle’s relative velocity Urei = v — g, where g 
is the velocity in the center of mass frame, is unchanged in magnitude (In'g/ = |nrei|) and 
with a randomly assigned angle such that all values of cos y, where y is the scattering angle, 
are equally likely. We assume, without loss of generality, that the direction of the particle’s 
velocity prior to the collision is aligned with the z axis. Note that we do not need to consider 
the electric field in this discussion. 

For a normal (“both”) collision, one can easily show [7] that scattering from a stationary 
neutral leads to an angle-dependent energy change given by 

2mmo 


ACboth — ^ 


-{cosx- 1 ), 


(m + mo)^ 

and a change of the ^-component of the momentum given by 

mo 


(Al) 


An 


2,both 


m + mo 


n(cos X — 1). 


(A2) 


Hence the average transfer of energy is 


('Ae)both — 2 


Acboth d cos X 


-1 


= 2mmoe/(m + mo)^ 


(A3) 

(A4) 


and the average transfer of momentum, for which the components perpendicular to 2 ; average 
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to zero, is 

(At;2)both = -moVz/{m + mo). (A5) 

For an “energy-only” collision, we select a random change in energy throngh a random 
“false” angle y, and calculate the energy change as if a normal collision through this angle 
had occurred. This preserves the average transfer of energy (Ae)energy = (Ae)both. We 
then discard the angle y, resetting the direction of motion of the particle to that before 
the collision. This has the unfortunate side effect of producing a change in momentum, 
contrary to the purpose of producing a collision with no momentum transfer. The change 
of momentum in this case is simply given by; 



and upon averaging over cos y we find; 

(An^)energy = -V 1 - ^ (l - (1 - 2a)^/^) 

m (3mo — m) 

3mo(m -|- mo) 

where a = 2mmo/(m -|- mo)^, and the second result is arrived at after some algebraic ma¬ 
nipulations, under the assumption that m < mo. 

For a “momentum only” collision, we perform a collision as normal using an angle y but set 
the post-collision energy to the pre-collision energy. Hence, it is obvious that (Ae)momentum=0 
but the modification of the energy also has an impact on the average momentum change, 
which we calculate below. Note that v' = vv', where v' is identical to that of a normal 
collision. Hence, 

V 

Koth ■ z ^-(mo cos y + m) (AT) 

m -|- mo 

implies that 
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From this, we find the change in the z-component of momentnm of a single collision 


An 


2,momentum 


V 


Z = viv 


1 ) 


(A9) 


and hence obtain an average momentum transfer of 


(Auz) 


momentum — 


cos X + — 

^ mo 




^l + 2^co.X+% 

^ 2 m\ 

3moJ 


— Id cos X 


(AlO) 


We must now combine these averages. In the Monte-Carlo simulation there are two 
distinct regimes depending on the value of r(e). If r(e) < 1, then the Uboth and cXenergy cross 
sections occur with frequencies such that; 

d{e) 

r{e)<i 

— ^spr(6) (A6)both “1“ i^sp(l r(6)) (A6)energy 
^sp(A6)both 

2mm, 


dt 


= 


{m + mo)^ 

where I'spiv) = nQvasp{v) would be the collision frequency in the absence of structure effects, 
and 


m 


d{v) 


dt 


r(d<i 


= mz/spr(e)(An)both + mu^p{3 - T{e)){Av) 


energy 


m + TUq 


mo 


(A12) 


where we have generalized the momentum transfer for any initial velocity direction, i.e. by 
replacing in equations (A5) and ( A10[ ) by v. In order to compare to the Boltzmann 
equations, we first note that <Jm{v) = <Jsp{v) for isotropic scattering and hence z^sp = 
and z/spr(e) = uA]"^. Hence, we can see that these rates of transfer are in agreement with 
the results for the Boltzmann equation (|^ and Q, excepting some factors of order m/mo 
that are neglected in our simulations and in the derivation of the Boltzmann equation. If 
r(e) > 1 then the “both” and “momentum only” cross sections occur with frequencies such 
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that; 


d{e) 


dt 


r{e)>i 

^sp('^*^)both “1“ ^sp(r(6) 1)(A6) momentum 

2mmo 


= z/. 


sp 


m + mo)^ 


(A13) 


and 


m 


d{v) 


dt 


rh)>i 


?7l^'sp('^^)both “1“ '^^sp(r(6) 1) momentum 

- --o-sp^.r(,) + o(-) 


. , . , (A14) 

m + mo mo 

and again we recover the Boltzmann eqnation results with negligible factors of m/mo. 

We note that it is possible to correct for these differences of the order of m/mo by 
modifying the cross sections cXenergy and cXmomentum- However, this would only be necessary 
if we considered systems in which the mass ratios were close to unity. 

Finally, we comment on the difference between the SSMC and WT methods. Wojcik and 
Tachiya [1] had mentioned that, in their particular calculation, they assumed the structure 
factor was mostly smaller than unity, which also implies that r(e) < 1. They did not deal 
with values greater than unity, instead choosing to set r(e) = min(r(e), 1). This means that 
the rates of energy and momentum transfer are unchanged when r(e) < 1 and continue to 


follow equations (All) and (A12). However, for r(e) > 1, this leads to an error in the WT 


method in the momentum transfer which is of the order: 

_mmon.pn 
m + mo 

There is no error in the energy transfer using the WT method. 


(A15) 


Appendix B: Benchmark Values for Transport Coefficients 
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the SSMC results. Quoted error is twice the standard error in the mean of the SSMC results. Comparisons between BE and SSMC 
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